Project 2

Code

using Plots;gr()

function RK4(dfunc,r,t,step)
    k1=dfunc(t,r)
    k2=dfunc(t+step/2,r.+k1*step/2)
    k3=dfunc(t+step/2,r.+k2*step/2)
    k4=dfunc(t+step,r.+k3*step)
    rn=r.+(k1.+2k2.+2k3.+k4)/6*step
    return rn
end
function run(func,r,t,step)
    xPoints=[]
    vPoints=[]
    for time in t
        xPoints=vcat(xPoints,r[1])
        vPoints=vcat(vPoints,r[2])
        r=RK4(func,r,time,step)
    end
    return xPoints,vPoints
end
function setRun(A,B,b,ω,stop,step=10^(-2))
    function func(t,r)
        x,v=r
        fX=v
        fV=B*cos(ω*t)-A*x^3-b*v
        return [fX,fV]
    end
    r=[3,0]
    ts=Array(range(0,stop,step=step))
    xs,vs=run(func,r,ts,step)
    return ts,xs,vs
end
function main()
    
    #For partA
    A=1;B=7;b=6;ω=1;stop=50pi
    #ts,xs,vs=setRun(A,B,b,ω,stop)
    #plot(xs,vs,label="v-x",xlabel="x",ylabel="v")
    #png(joinpath(@__DIR__,"Project/partA.png"))

    #For partB
    #for (B,b) in [(7,6),(7,0.6),(10,0.05),(7,0.01)]
    #    ts,xs,vs=setRun(A,B,b,ω,stop)
    #    plot(ts,xs,label="position",xlabel="t",ylabel="x")
    #    png(joinpath(@__DIR__,"Project/partB/position_$B&$b.png"))
    #    plot(xs,vs,label="v-x",xlabel="x",ylabel="v")
    #    png(joinpath(@__DIR__,"Project/partB/phase_$B&$b.png"))
    #end

    #For partC
    N=1000
    b=0.01;step=pi/180;stop=360*step*N
    ts,xs,vs=setRun(A,B,b,ω,stop,step)
    if mod(length(vs),2)!=0
        pop!.([vs,xs])
    end
    vs,xs=reshape.([vs,xs],360,:)
    anim=@animate for off=1:360
        V=vs[off,:]
        X=xs[off,:]
        scatter(X,V,xlims=[-4,4],ylims=[-6,6],xlabel="x",ylabel="v",legend=false,title="phase offset: $off °",markersize=0.8,color=:black)
    end
    gif(anim,joinpath(@__DIR__,"Project/partC.gif"),fps=10)
end
main()

运行结果及分析

Part A

Figure 1 Part A

由图中我们可以清楚看出,该运动最后归于稳定解.

Part B

Figure 2 相图(B=7,b=6)

Figure 3 位图(B=7,b=6)

Figure 4 相图(B=7,b=0.6)

Figure 5 位图(B=7,b=0.6)

Figure 6 相图(B=7,b=0.01)

Figure 7 位图(B=7,b=0.01)

Figure 8 相图(B=10,b=0.05)

Figure 9 位图(B=10,b=0.05)

我们可以看出除了(B,b)=(7,0.01)的组合为混沌情况之外,其余的组合均为稳定解.
而(B,b)=(10,0.05)为分叉.

Part C

修改推导的时间间隔使得每个间隔相位变换为1°,同时设置总推导周期数为1000,每个庞加莱截面的点数为1000,最后所得结果为

Figure 10 庞加莱截面